Annual Sea Level Rankings

In this notebook, we’ll be creating a table and a time series plot of sea level rankings at Malakal from 1993-2024. Looking at extreme sea levels by ranking our top events helps us quickly put high and low water levels into historical context with relatively little coding effort.

Setup

As with previous sections, we first need to import the necessary libraries, and establish our input/output directories. Note that this notebook is also largely a repeat of the anomaly notebook.

# import necessary libraries
import numpy as np
import xarray as xr
import datetime as dt
from pathlib import Path
import pandas as pd
import os
import os.path as op
import sys

# data retrieval libraries
import requests
from urllib.request import urlretrieve #used for downloading files
import json
import copernicusmarine

# data processing libraries
from scipy import stats
import seaborn as sns
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

from myst_nb import glue #used for figure numbering when exporting to LaTeX

sys.path.append("../../../functions")
from data_downloaders import download_oni_index
from sea_level_func import detect_enso_events

Now we’ll set up our input and output directories.

data_dir = Path('../../../data')
path_figs = "../../../matrix_cc/figures"
data_dir = Path(data_dir,'sea_level')
output_dir = data_dir 

# Create the output directory if it does not exist
output_dir.mkdir(exist_ok=True)
data_dir.mkdir(exist_ok=True)

# Also let's just make our figure sizes the same:
plt.rcParams['figure.figsize'] = [10, 4]  # Set a default figure size for the notebook

Retrieve the Tide Station(s) Data Set(s)

Next, we’ll access data from the UHSLC. The Malakala tide gauge has the UHSLC ID: 7. We will import the netcdf file into our current data directory, and take a peek at the dataset. We will also import the datums for this location.

uhslc_id = 7
# download the hourly data
fname = f'h{uhslc_id:03}.nc' # h for hourly, d for daily

url = "https://uhslc.soest.hawaii.edu/data/netcdf/fast/hourly/" 

path = os.path.join(data_dir, fname)

if not os.path.exists(path):
    urlretrieve(os.path.join(url, fname), path) 

rsl = xr.open_dataset(data_dir / fname)

# remove the trailing zero from record_id. This zero is added to the record_id to make it unique if the station has multiple entries
rsl['record_id'] =(rsl['record_id']/10).astype(int)

# change station_name to string
for col in ['station_name', 'station_country']:
    rsl[col] = rsl[col].astype(str)
rsl

<xarray.Dataset> Size: 6MB
Dimensions:               (record_id: 1, time: 486849)
Coordinates:
  * time                  (time) datetime64[ns] 4MB 1969-05-18T15:00:00 ... 2...
  * record_id             (record_id) int64 8B 7
Data variables:
    sea_level             (record_id, time) float32 2MB ...
    lat                   (record_id) float32 4B ...
    lon                   (record_id) float32 4B ...
    station_name          (record_id) <U7 28B 'Malakal'
    station_country       (record_id) <U5 20B 'Palau'
    station_country_code  (record_id) float32 4B ...
    uhslc_id              (record_id) int16 2B ...
    gloss_id              (record_id) float32 4B ...
    ssc_id                (record_id) |S4 4B ...
    last_rq_date          (record_id) datetime64[ns] 8B ...
Attributes:
    title:                  UHSLC Fast Delivery Tide Gauge Data (hourly)
    ncei_template_version:  NCEI_NetCDF_TimeSeries_Orthogonal_Template_v2.0
    featureType:            timeSeries
    Conventions:            CF-1.6, ACDD-1.3
    date_created:           2025-01-27T14:35:24Z
    publisher_name:         University of Hawaii Sea Level Center (UHSLC)
    publisher_email:        philiprt@hawaii.edu, markm@soest.hawaii.edu
    publisher_url:          http://uhslc.soest.hawaii.edu
    summary:                The UHSLC assembles and distributes the Fast Deli...
    processing_level:       Fast Delivery (FD) data undergo a level 1 quality...
    acknowledgment:         The UHSLC Fast Delivery database is supported by ...

xarray.Dataset

• Dimensions:
  • record_id: 1
  • time: 486849
• Coordinates: (2)
  • time
    (time)
    datetime64[ns]
    1969-05-18T15:00:00 ... 2024-11-...

    long_name :

    time

    axis :

    T

    array(['1969-05-18T15:00:00.000000000', '1969-05-18T16:00:00.028800000',
           '1969-05-18T16:59:59.971200000', ..., '2024-11-30T21:00:00.000000000',
           '2024-11-30T22:00:00.028800000', '2024-11-30T22:59:59.971200000'],
          shape=(486849,), dtype='datetime64[ns]')

  • record_id
    (record_id)
    int64
    7

    array([7])

• Data variables: (10)
  • sea_level
    (record_id, time)
    float32
    ...

    long_name :

    relative sea level

    units :

    millimeters

    source :

    in situ tide gauge water level observations

    platform :

    station_name, station_country, station_country_code, uhslc_id, gloss_id, ssc_id

    [486849 values with dtype=float32]

  • lat
    (record_id)
    float32
    ...

    standard_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    valid_min :

    -90

    valid_max :

    90

    [1 values with dtype=float32]

  • lon
    (record_id)
    float32
    ...

    standard_name :

    longitude

    units :

    degrees_east

    axis :

    X

    valid_min :

    0

    valid_max :

    360

    [1 values with dtype=float32]

  • station_name
    (record_id)
    <U7
    'Malakal'

    long_name :

    station name

    array(['Malakal'], dtype='<U7')

  • station_country
    (record_id)
    <U5
    'Palau'

    long_name :

    station country (ISO 3166-1)

    array(['Palau'], dtype='<U5')

  • station_country_code
    (record_id)
    float32
    ...

    long_name :

    station country code (ISO 3166-1 numeric)

    comment :

    These are 3-digit country codes (e.g., 003) stored as integers.

    [1 values with dtype=float32]

  • uhslc_id
    (record_id)
    int16
    ...

    long_name :

    unique station ID number used by the University of Hawaii Sea Level Center (UHSLC)

    [1 values with dtype=int16]

  • gloss_id
    (record_id)
    float32
    ...

    long_name :

    unique station ID number used by the WMO/IOC Global Sea Level Observing System (GLOSS)

    [1 values with dtype=float32]

  • ssc_id
    (record_id)
    |S4
    ...

    long_name :

    unique station ID code in the Sealevel Station Catalog (SSC) produced by the WMO/IOC Sea Level Monitoring Facility (VLIZ)

    comment :

    Note that SSC IDs vary in length. The IDs are padded with space characters to produce the retangular character matrix provided here.

    [1 values with dtype=|S4]

  • last_rq_date
    (record_id)
    datetime64[ns]
    ...

    long_name :

    date of last Research Quality sea level value

    comment :

    Values in last_rq_date correspond to an identical value in the time vector. Dates less than or equal to last_rq_date for a given station correspond to Research Quality (RQ) data. Dates after last_rq_date correspond to Fast Delivery (FD) data. Missing values indicate no RQ data for a station.

    [1 values with dtype=datetime64[ns]]

• Indexes: (2)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex([       '1969-05-18 15:00:00', '1969-05-18 16:00:00.028800',
                   '1969-05-18 16:59:59.971200',        '1969-05-18 18:00:00',
                   '1969-05-18 19:00:00.028800', '1969-05-18 19:59:59.971200',
                          '1969-05-18 21:00:00', '1969-05-18 22:00:00.028800',
                   '1969-05-18 22:59:59.971200',        '1969-05-19 00:00:00',
                   ...
                   '2024-11-30 13:59:59.971200',        '2024-11-30 15:00:00',
                   '2024-11-30 16:00:00.028800', '2024-11-30 16:59:59.971200',
                          '2024-11-30 18:00:00', '2024-11-30 19:00:00.028800',
                   '2024-11-30 19:59:59.971200',        '2024-11-30 21:00:00',
                   '2024-11-30 22:00:00.028800', '2024-11-30 22:59:59.971200'],
                  dtype='datetime64[ns]', name='time', length=486849, freq=None))

  • record_id
    PandasIndex

    PandasIndex(Index([7], dtype='int64', name='record_id'))

• Attributes: (11)

  title :

  UHSLC Fast Delivery Tide Gauge Data (hourly)

  ncei_template_version :

  NCEI_NetCDF_TimeSeries_Orthogonal_Template_v2.0

  featureType :

  timeSeries

  Conventions :

  CF-1.6, ACDD-1.3

  date_created :

  2025-01-27T14:35:24Z

  publisher_name :

  University of Hawaii Sea Level Center (UHSLC)

  publisher_email :

  philiprt@hawaii.edu, markm@soest.hawaii.edu

  publisher_url :

  http://uhslc.soest.hawaii.edu

  summary :

  The UHSLC assembles and distributes the Fast Delivery (FD) dataset of hourly- and daily-averaged tide gauge water-level observations. Tide gauge operators, or data creators, provide FD data to UHSLC after a level 1 quality assessment (see processing_level attribute). The UHSLC provides an independent quality assessment of the time series and makes FD data available within 4-6 weeks of collection. This is a "fast" turnaround time compared to Research Quality (RQ) data, which are available on an annual cycle after a level 2 quality assessment. RQ data replace FD data in the data stream as the former becomes available. This file contains hybrid time series composed of RQ data when available with FD data appended to the end of each RQ series.

  processing_level :

  Fast Delivery (FD) data undergo a level 1 quality assessment (e.g., unit and timing evaluation, outlier detection, combination of multiple channels into a primary channel, etc.). In this file, FD data are appended to Research Quality (RQ) data that have received a level 2 quality assessment (e.g., tide gauge datum evaluation, assessment of level ties to tide gauge benchmarks, comparison with nearby stations, etc.).

  acknowledgment :

  The UHSLC Fast Delivery database is supported by the National Oceanic and Atmospheric Administration (NOAA) Office of Climate Observations (OCO).

Adjust to MHHW

def get_MHHW_uhslc_datums(id, datumname): 
    
    url = 'https://uhslc.soest.hawaii.edu/stations/TIDES_DATUMS/fd/LST/fd'+f'{int(id):03}'+'/datumTable_'+f'{int(id):03}'+'_mm_GMT.csv'
    datumtable = pd.read_csv(url)
    datum = datumtable[datumtable['Name'] == datumname]['Value'].values[0]
    # ensure datum is a float
    datum = float(datum)
    return datum, datumtable

# Adjust to MHHW in meters and update attributes
mhhw = get_MHHW_uhslc_datums(uhslc_id, 'MHHW')[0]
sea_level_attrs = rsl['sea_level'].attrs.copy()
rsl['sea_level'] = (rsl['sea_level'] - mhhw) / 1000
rsl['sea_level'].attrs = sea_level_attrs  # restore all previous attributes
rsl['sea_level'].attrs['datum'] = 'MHHW (1983-2001)'  # update/add as needed
rsl['sea_level'].attrs['units'] = 'meters'

Resample to monthly

For background plotting purposes. Here we’ll extract the mean, min and max of the monthly MHHW sea level.

# get min and max for each month by resampling
rsl_monthly_min = rsl.resample(time='1ME').min()
rsl_monthly_max = rsl.resample(time='1ME').max()
rsl_monthly_mean = rsl.resample(time='1ME').mean()

Define functions

The first will get us the top 10 high and low hourly sea level values for a given station. To ensure that we are tracking unique high and low water events, we ensure that the hourly maxima and minima are separated by at least 3 days.

def get_top_ten(rsl, record_id, mode='max'):
    # Convert data to a pandas Series
    sea_level_series = rsl.sea_level.sel(record_id=record_id).to_series()

   # Select top 100 values based on the mode
    if mode == 'max':
        top_values = sea_level_series.nlargest(100)
    elif mode == 'min':
        top_values = sea_level_series.nsmallest(100)
    else:
        raise ValueError('mode must be either "max" or "min"')

    # Filter to find unique events spaced by at least 3 days
    filtered_dates = []
    top_10_values = pd.Series()

    for date, value in top_values.items():
        if all(abs((date - pd.to_datetime(added_date)).days) > 3 for added_date in filtered_dates):
            filtered_dates.append(date)
            top_10_values[date] = value
        if len(filtered_dates) == 10:
            break
    rank = np.arange(1,11)

    station_name = str(rsl['station_name'].sel(record_id=record_id).values)
    # station_name = str(rsl['station_name'].isel(record_id=record_id).values)
    # record_id = rsl['record_id'].isel(record_id=record_id).values

    top_10_values = pd.DataFrame({'rank': rank, 'date': top_10_values.index, 'sea level (m)': top_10_values.values})  
    top_10_values['station_name'] = station_name
    top_10_values['record_id'] = record_id
    top_10_values['type'] = mode

    #round the date to the nearest hour
    top_10_values['date'] = top_10_values['date'].dt.round('h')

    return top_10_values

Make a Table

We’ll use the function we defined above to make a table of the top 10 ranked high and low water events. We’ll add the ONI values to this table, which will come in handy later.

def get_top_10_table(rsl,record_id):
    # make a table of the top 10 values, sorted by size and with date
    top_10_values_max = get_top_ten(rsl, record_id, mode='max')
    top_10_values_min = get_top_ten(rsl, record_id, mode='min')

    top_10_table = pd.concat([top_10_values_max,top_10_values_min])

    # cross reference the dates with the oni data to see if they are during an El Nino or La Nina event
    p_data = 'https://psl.noaa.gov/data/correlation/oni.data'
    oni = download_oni_index(p_data)
    oni = detect_enso_events(oni)

    #drop the La Nina and El Nino columns
    oni = oni.drop(columns=['La Nina', 'El Nino'])

    #Extract ONI values for the corresponding dates of top_10_table
    oni_val = oni.reindex(top_10_table['date'], method='nearest')
    
    # add the oni values to the top_10_table
    top_10_table = pd.merge(top_10_table, oni_val, left_on='date', right_index=True)

    return top_10_table

Note that because we started this analysis with hourly data, we may not be fully aligned with max/min statistics from other sources. Hourly data will dampen both the minima and maxima due to averaging. Note also that there are other extreme water level rankings that may be offset by a different MHHW datum. Here, we used the 1983-2001 NTDE, corresponding to the datum used most often by NOAA CO-OPS.

For other ranking tables, check out the NASA Flooding Analysis Tool or the Sea Level Calculator

station_name = str(rsl['station_name'].sel(record_id=uhslc_id).values)

top_10_table = get_top_10_table(rsl, uhslc_id)

station_name

/Users/laurac/Library/Mobile Documents/com~apple~CloudDocs/Projects/CC_indicators/CC_indicators/notebooks_historical/ocean/1_sea_level/../../../functions/data_downloaders.py:242: FutureWarning: The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead
  oni = pd.read_csv(

'Malakal'

Display by rank

Now we’ll start to make the table a little more presentable by labeling the columns properly, and organizing events by rank.

# subdivide the data into two columns for max and min
# Filter out the 'max' and 'min' types into separate DataFrames
max_top_10 = top_10_table[top_10_table['type'] == 'max'].reset_index(drop=True)
min_top_10 = top_10_table[top_10_table['type'] == 'min'].reset_index(drop=True)

# make table of highest and lowest values by rank, with columns: rank, highest sea level (m), date, lowest sea level (m), date
top_10_display = pd.DataFrame({'Rank':max_top_10['rank'],
                             'Highest':max_top_10['sea level (m)'].apply(lambda x: float(f"{x:.3g}")),
                             'Highest Date':max_top_10['date'].dt.strftime('%Y-%m-%d %H:%M'),
                             'Highest ONI Mode':max_top_10['ONI Mode'],
                             'Lowest':min_top_10['sea level (m)'].apply(lambda x: float(f"{x:.3g}")),
                             'Lowest Date':min_top_10['date'].dt.strftime('%Y-%m-%d %H:%M'),
                             'Lowest ONI Mode':min_top_10['ONI Mode'],
                             'Zone': 'GMT',
                             'ONI max':max_top_10['ONI'],
                             'ONI min':min_top_10['ONI']})

Style the table

Now we’ll add some style. The table will be colored by the ONI, if the event is classified as an El Nino or La Nina.

top_10_display

	Rank	Highest	Highest Date	Highest ONI Mode	Lowest	Lowest Date	Lowest ONI Mode	Zone	ONI max	ONI min
0	1	0.625	2013-06-23 22:00	Neutral	-2.19	1992-01-19 17:00	El Nino	GMT	-0.40	1.63
1	2	0.613	2016-10-16 10:00	La Nina	-2.18	1986-12-31 17:00	El Nino	GMT	-0.69	1.23
2	3	0.589	2016-05-07 22:00	Neutral	-2.17	1982-12-31 18:00	El Nino	GMT	0.39	2.18
3	4	0.582	2011-09-29 11:00	La Nina	-2.12	2002-12-04 17:00	El Nino	GMT	-1.01	1.14
4	5	0.581	1998-10-07 11:00	La Nina	-2.11	1995-01-01 17:00	El Nino	GMT	-1.35	0.96
5	6	0.576	2011-03-21 23:00	La Nina	-2.11	1970-01-08 18:00	El Nino	GMT	-0.73	0.51
6	7	0.575	2024-10-19 11:00	Neutral	-2.10	2019-01-21 17:00	El Nino	GMT	-0.37	0.72
7	8	0.572	2021-03-29 23:00	La Nina	-2.10	1983-01-28 17:00	El Nino	GMT	-0.66	1.92
8	9	0.564	1989-01-10 12:00	La Nina	-2.10	1982-12-02 18:00	El Nino	GMT	-1.69	2.23
9	10	0.561	2007-09-26 22:00	La Nina	-2.10	1987-01-29 17:00	El Nino	GMT	-1.34	1.19

Show code cell source Hide code cell source

# # Define a function that returns styles for both 'Highest' and 'Lowest' based on respective ONI Modes
# # Normalize the ONI values
import matplotlib.colors as mcolors
norm = mcolors.Normalize(vmin=-3, vmax=3)

# Define the color maps
cmap = sns.color_palette('coolwarm', as_cmap=True)

def style_oni_based(row):
    styles = {}
    # Color for 'Highest' based on 'Highest ONI Mode'
    if row['Highest ONI Mode'] == 'El Nino':
        styles['Highest'] = f'background-color: {mcolors.rgb2hex(cmap(norm(row["ONI max"])))}'
    elif row['Highest ONI Mode'] == 'La Nina':
        styles['Highest'] = f'background-color: {mcolors.rgb2hex(cmap(norm(row["ONI max"])))}'
    else:
        styles['Highest'] = 'color: black'  # Neutral or undefined
    # Color for 'Lowest' based on 'Lowest ONI Mode'
    if row['Lowest ONI Mode'] == 'El Nino':
        styles['Lowest'] = f'background-color: {mcolors.rgb2hex(cmap(norm(row["ONI min"])))}'
    elif row['Lowest ONI Mode'] == 'La Nina':
        styles['Lowest'] = f'background-color: {mcolors.rgb2hex(cmap(norm(row["ONI min"])))}'
    else:
        styles['Lowest'] = 'color: black'  # Neutral or undefined
    return pd.Series(styles)    

# Define the table style
styles_col = [
    {'selector': 'th', 'props': [('border-bottom', '2px solid #000'), ('font-weight', 'bold'), ('text-align', 'center')]},
    {'selector': '.dataframe th:nth-child(3)', 'props': [('padding-right', '100px')]}]

# Apply the styling function to the DataFrame
top_10_d = top_10_display.style.\
                            apply(style_oni_based, axis=1).\
                            hide(subset=['Highest ONI Mode', 'Lowest ONI Mode','ONI max','ONI min'], axis=1).\
                            hide().\
                            set_table_attributes(styles_col)


# Display the styled DataFrame



tablename = 'SL_top_10_table_' + station_name + '.csv'
# save table to output folder
top_10_display.to_csv(output_dir / tablename, index=False)

top_10_d

Rank	Highest	Highest Date	Lowest	Lowest Date	Zone
1	0.625000	2013-06-23 22:00	-2.190000	1992-01-19 17:00	GMT
2	0.613000	2016-10-16 10:00	-2.180000	1986-12-31 17:00	GMT
3	0.589000	2016-05-07 22:00	-2.170000	1982-12-31 18:00	GMT
4	0.582000	2011-09-29 11:00	-2.120000	2002-12-04 17:00	GMT
5	0.581000	1998-10-07 11:00	-2.110000	1995-01-01 17:00	GMT
6	0.576000	2011-03-21 23:00	-2.110000	1970-01-08 18:00	GMT
7	0.575000	2024-10-19 11:00	-2.100000	2019-01-21 17:00	GMT
8	0.572000	2021-03-29 23:00	-2.100000	1983-01-28 17:00	GMT
9	0.564000	1989-01-10 12:00	-2.100000	1982-12-02 18:00	GMT
10	0.561000	2007-09-26 22:00	-2.100000	1987-01-29 17:00	GMT

#make this a pretty table with great_tables
from great_tables import GT, md, html

# replace all spaces with underscores in top_10_display column names
top_10_display.columns = top_10_display.columns.str.replace(' ', '_')

outFilename = 'SL_top_10_table_' + station_name + '.png'
outFilepath = output_dir / outFilename

# Convert the DataFrame to a great_table
(
    GT(top_10_display)
    .tab_header(title='Top 10 Sea Level Events', subtitle='Station: '+station_name)
    .tab_spanner(label='HIGHEST', columns=['Highest', 'Highest_Date', 'Highest_ONI_Mode'])
    .tab_spanner(label='LOWEST', columns=['Lowest', 'Lowest_Date', 'Lowest_ONI_Mode'])
    .cols_hide(columns=['ONI_max','ONI_min','Zone'])
    .cols_align('center')
    .cols_label(Highest=html('Water Level <br>(m, MHHW)'), Highest_Date='Date', Highest_ONI_Mode='ONI Mode',
                Lowest=html('Water Level <br>(m, MHHW)'), Lowest_Date='Date', Lowest_ONI_Mode='ONI Mode')
                .tab_source_note(source_note='Source: UHSLC Hourly Water Level, https://psl.noaa.gov/data/correlation/oni.data')
    .save(outFilepath)
)

Top 10 Sea Level Events
Station: Malakal
Rank	HIGHEST			LOWEST
	Water Level (m, MHHW)	Date	ONI Mode	Water Level (m, MHHW)	Date	ONI Mode
1	0.625	2013-06-23 22:00	Neutral	-2.19	1992-01-19 17:00	El Nino
2	0.613	2016-10-16 10:00	La Nina	-2.18	1986-12-31 17:00	El Nino
3	0.589	2016-05-07 22:00	Neutral	-2.17	1982-12-31 18:00	El Nino
4	0.582	2011-09-29 11:00	La Nina	-2.12	2002-12-04 17:00	El Nino
5	0.581	1998-10-07 11:00	La Nina	-2.11	1995-01-01 17:00	El Nino
6	0.576	2011-03-21 23:00	La Nina	-2.11	1970-01-08 18:00	El Nino
7	0.575	2024-10-19 11:00	Neutral	-2.1	2019-01-21 17:00	El Nino
8	0.572	2021-03-29 23:00	La Nina	-2.1	1983-01-28 17:00	El Nino
9	0.564	1989-01-10 12:00	La Nina	-2.1	1982-12-02 18:00	El Nino
10	0.561	2007-09-26 22:00	La Nina	-2.1	1987-01-29 17:00	El Nino
Source: UHSLC Hourly Water Level, https://psl.noaa.gov/data/correlation/oni.data

Make Timeseries Plots

Now we’ll plot this data up in the time domain.

Interactive Plots

This version is interactive. Bad for making pdf reports though.

import plotly.graph_objects as go

def make_plotlyFigure(rsl_monthly_mean, rsl_monthly_max, rsl_monthly_min, top_10_table, rsl_subset, record_id):
    figly = go.Figure()
    top_10_table = get_top_10_table(rsl_subset, record_id)
    
    # Select data for the current station
    max_sl = rsl_monthly_max['sea_level'].sel(record_id=record_id)
    min_sl = rsl_monthly_min['sea_level'].sel(record_id=record_id)
    mean_sl = rsl_monthly_mean['sea_level'].sel(record_id=record_id)

    x = mean_sl.time - np.timedelta64(15, 'D')  # add 15 days to make the line in the middle of the month


    # Drop NaN values and get the corresponding time coordinates
    max_sl_clean = max_sl.dropna(dim='time')
    min_sl_clean = min_sl.dropna(dim='time')
    mean_sl_clean = mean_sl.dropna(dim='time')


    # Add traces
    figly.add_trace(go.Scatter(x=x, y=max_sl_clean, 
                    mode='lines', line_color='rgba(0,176,246,0.2)', name='Monthly Max'))

    figly.add_trace(go.Scatter(x=x, y=min_sl_clean, 
                    mode='lines', line_color='rgba(0,176,246,0.2)', fill='tonexty', 
                    fillcolor='rgba(0,176,246,0.2)', name='Range (Monthly Max/Min)'))

    figly.add_trace(go.Scatter(x=x, y=mean_sl_clean, 
                    mode='lines', line_color='rgba(0,176,246,1)', name='Monthly Mean'))

    # Add El Niño markers
    figly.add_trace(go.Scatter(x=top_10_table[top_10_table['ONI Mode']=='El Nino'].date,
                                y=top_10_table[top_10_table['ONI Mode']=='El Nino']['sea level (m)'],
                                mode='markers', marker=dict(color='red', size=12, symbol='star'), name='El Niño',
                                hoverinfo='none'))
    
    # Add La Niña markers
    figly.add_trace(go.Scatter(x=top_10_table[top_10_table['ONI Mode']=='La Nina'].date,
                                y=top_10_table[top_10_table['ONI Mode']=='La Nina']['sea level (m)'],
                                mode='markers', marker=dict(color='blue', size=12, symbol='circle'), name='La Niña',
                                hoverinfo='none'))
    
    # Add highest/lowest events markers
    figly.add_trace(go.Scatter(x=top_10_table.date, y=top_10_table['sea level (m)'], 
                                mode='markers', marker=dict(color='orange', size=6, symbol='circle'), name='Highest/Lowest Events',
                                hovertemplate='%{x}, %{y:.2f} m<extra></extra>'))
    
    # Update layout
    figly.update_layout(
        title=str(station_name),
        xaxis_title='Time',
        yaxis_title=f"{rsl_subset['sea_level'].attrs['long_name']} [{rsl_subset['sea_level'].attrs['units']}]",
        showlegend=True,
        # width=500,  # Set figure width
        # height=400,  # Set figure height
        autosize=True,  # Turn off autosize
        plot_bgcolor='white',
        xaxis=dict(showgrid=True, gridwidth=0.5, gridcolor='rgb(200,200,200)', range=['1993-01-01', '2024-12-31']),
        yaxis=dict(showgrid=True, gridwidth=0.5, gridcolor='rgb(200,200,200)'),
        legend=dict(
            orientation='h',  # Horizontal orientation
            y=-0.3,  # Position below the plot
            x=0.5,  # Centered horizontally
            xanchor='center',
            yanchor='top',
            font=dict(size=10)  # Smaller font size for legend items
        ),
        margin=dict(l=50, r=50, t=50, b=100),  # Increase bottom margin to make room for the legend
        title_x=0.5,  # Center the title
        title_font_size=20,
        title_font_family='Arial',
        title_font_color='black'
    )


    
    # figly.show()
    return figly

A note: the plotly plot will not render in the jupyter book build, because the original data to build it is stored elsewhere (not on github). If you are replicating this notebook, it should work fine for you. In the meantime, the pre-built html is available.

figly = make_plotlyFigure(rsl_monthly_mean, 
        rsl_monthly_max, rsl_monthly_min, 
        top_10_table, rsl,uhslc_id)

# change size of the plot
# figly.update_layout(autosize=True)
figly.update_layout(width=800, height=500,autosize=False)
# save the figure to the output folder
figname = 'SL_top_10_plot_' + station_name + '.html'

# remove commas from figname
figname = figname.replace(', ','_')


figly

/Users/laurac/Library/Mobile Documents/com~apple~CloudDocs/Projects/CC_indicators/CC_indicators/notebooks_historical/ocean/1_sea_level/../../../functions/data_downloaders.py:242: FutureWarning:

The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead

Static Plots

Plotting them all so we can decide what to do with them. Click “show code cell output” to see all stations.

Show code cell source Hide code cell source

# static version
top_10_table = get_top_10_table(rsl,uhslc_id)
# get top_10_table for each record
fig,ax = plt.subplots(1,1,figsize=(10,3))

x = rsl_monthly_mean.time - np.timedelta64(15,'D') # add 15 days to make the line in the middle of the month

# station_name = str(rsl_subset['station_name'].sel(station_id=station_id).values)
ax.set_title(station_name)
ax.set_ylabel(f'{rsl["sea_level"].attrs["long_name"]} [{rsl["sea_level"].attrs["units"]}]')

# Plot the max, min, and mean lines
# Plot the max, min, and mean lines
ax.plot(x, rsl_monthly_max['sea_level'].sel(record_id=uhslc_id), color=(0/255, 176/255, 246/255, 0.2))
ax.plot(x, rsl_monthly_min['sea_level'].sel(record_id=uhslc_id), color=(0/255, 176/255, 246/255, 0.2))
ax.fill_between(x, rsl_monthly_max['sea_level'].sel(record_id=uhslc_id), 
                 rsl_monthly_min['sea_level'].sel(record_id=uhslc_id), color=(0/255, 176/255, 246/255, 0.2), label='Range (Monthly Max/Min)')
ax.plot(x, rsl_monthly_mean['sea_level'].sel(record_id=uhslc_id), color=(0/255, 176/255, 246/255, 1), label='Monthly Mean')
# Plot the markers for La Niña, El Niño, and other events
ax.scatter(top_10_table[top_10_table['ONI Mode']=='La Nina'].date, top_10_table[top_10_table['ONI Mode']=='La Nina']['sea level (m)'], color='blue', s=120, label='La Niña', marker='o')
ax.scatter(top_10_table[top_10_table['ONI Mode']=='El Nino'].date, top_10_table[top_10_table['ONI Mode']=='El Nino']['sea level (m)'], color='red', s=120, label='El Niño', marker='*')
ax.scatter(top_10_table.date, top_10_table['sea level (m)'], color='orange', s=20, label='Highest/Lowest Events', marker='o')
    
# Setting the x-axis limits
ax.set_xlim(pd.Timestamp('1993-01-01'), pd.Timestamp('2024-12-31'))

# Improve the appearance
ax.grid(True, color='gray', alpha=0.5 ,linewidth=0.5)
ax.legend(loc='upper center', bbox_to_anchor=(0.5, -0.2), ncol=3, frameon=False, fontsize=10)


# Save the plot to a file
figname = 'SL_rankings_' + station_name.replace(' ', '_') + '.png'
fig.savefig(output_dir / figname, dpi=300, bbox_inches='tight')


glue("SL_rankings",fig,display=False)

Show code cell output Hide code cell output

/Users/laurac/Library/Mobile Documents/com~apple~CloudDocs/Projects/CC_indicators/CC_indicators/notebooks_historical/ocean/1_sea_level/../../../functions/data_downloaders.py:242: FutureWarning:

The 'delim_whitespace' keyword in pd.read_csv is deprecated and will be removed in a future version. Use ``sep='\s+'`` instead

Fig. 11 This is a static version of the previous figure, which can be used for print (or static websites).

Citations